COMMENT ⊗   VALID 00015 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002	LINE - "hardware"-type routines for edges, lines, vertices.
C00005 00003	_ SUMS
C00009 00004	_ MALIN, ANGDIR, ANGLE
C00012 00005	_ ANGLEN, LINFIT
C00015 00006	_ LINFIT cont
C00017 00007	_ LINFIT cont, DNEW, ICON
C00020 00008	_ SORTED
C00022 00009	_ SORTED cont.
C00024 00010	_ SORTED cont.
C00026 00011	_ SORTED cont.
C00029 00012	_ WEIGHV
C00032 00013	_ MALI, PLDIS
C00034 00014	_ LDIST, SKAR1, LEKV
C00036 00015	_ REKOP, INREK
C00039 ENDMK
C⊗;
COMMENT LINE - "hardware"-type routines for edges, lines, vertices.;

ENTRY SUMS,MALIN,ANGLEN,LINFIT,DNEW,ICON,SORTED,MALI,WEIGHV,PLDIS,
	LDIST,SKAR1,LEKV,REKOP,ANGDIR,INREK;

BEGIN "LINE"

DEFINE QI="INTEGER",
	QR="REAL",
	QRI="REFERENCE INTEGER",
	QRR="REFERENCE REAL",
	QEP="EXTERNAL SIMPLE PROCEDURE",
	QEIP="EXTERNAL SIMPLE INTEGER PROCEDURE",
	QERP="EXTERNAL SIMPLE REAL PROCEDURE",
	QFOP="FORWARD INTERNAL SIMPLE PROCEDURE",
	QFOIP="FORWARD INTERNAL SIMPLE INTEGER PROCEDURE",
	QFORP="FORWARD INTERNAL SIMPLE REAL PROCEDURE",
	_="COMMENT",
	LOOP(I,J,K,L)="FOR I←J STEP L UNTIL K DO",
	SAFEX="SAFE";

REAL A1;
SAFEX INTERNAL INTEGER ARRAY VERLIN[1:20];
EXTERNAL INTEGER NOEPA,NOL,NOV,IFREEL,NOBAL,LDATE,ILLL,ILFL;
EXTERNAL REAL RDEP,RMEDA,SHRINK,RDDP,RDNP,RMSD,RMLG;
INTERNAL REAL A11, A12, A21, A22, X00, Y00;
SAFEX EXTERNAL INTEGER ARRAY LE[0:1],LEDG1,LEDG2,LCREDE[1:1];
SAFEX EXTERNAL REAL ARRAY EAX,EAY,EBX,EBY[0:1],XLCOR,YLCOR,CXL,CYL,
	CCL,RLEN,ANGARG[1:1];

	QEP TELL(STRING S);
	QEP UNTELL;
	QEIP KARN(QR X1,Y1,X2,Y2,X3,Y3,X4,Y4; QRR X,Y; QRI IX1,IX2,IP1,IP2;
		  QRR R1,R2; QI IC; QR WI);
	QEP KOT(QR X1,Y1,X2,Y2,X3,Y3);
	QFOP ANGLEN(QI LL);
	QFOIP ICON(QI I,J,K);
	QFORP  DNEW(QR A,B,C,D,E,F,G);
	QEIP  LVNEXT(QI I,J);
	QEP RETLIN(QI I);
	QEIP MSCVCO(QI I,J,K);
	QFOP LEKV(QR X1,Y1,X2,Y2; QRR A,B,C);
	QERP SQRT(QR R);
	QERP COS(QR R);
	QERP ATAN2(QR R,S);
	QERP SIGN(QR R,S);
	QERP AMOD(QR R,S);
_ SUMS;

_	Adds one or more edge-pairs in the various sums
	for the least square fit.
	Returns new mean square deviation, if tolerable,
	else -100.;

INTERNAL SIMPLE REAL PROCEDURE SUMS(INTEGER ILO,IHI1;
		REFERENCE REAL CX,CY,CC; REAL TOL);
	BEGIN "SUMS"
	LABEL L10,L100,L101,L1,L11,L40,L51,L4,L50,L52;
	INTEGER IHI,IHI2;
	REAL RRET,E11,E12,E21,E22,SX,SY,SX2,SY2,SXY,CN,F,G,S1,S2,S3,S4,ROT,
		SQ1,SQ2,SQT;
	RRET←-100.;
	IF IHI1>0 THEN GO L10;
	IHI←ILO;
	IHI2←-IHI1;
	GO L100;

L10:	IHI←IHI1;
L100:	E11←EAX[IHI];
	E21←EBX[IHI];
	E12←EAY[IHI];
	E22←EBY[IHI];
	IF IHI≠ILO THEN GO L1;
	SX←E11+E21;
	SY←E12+E22;
	SX2←E11*E11+E21*E21;
	SY2←E12*E12+E22*E22;
	SXY←E11*E12+E21*E22;
	LEKV(E11,E12,E21,E22,CX,CY,CC);
	IF IHI1<0 THEN GO L101;
	RETURN(0.);

L101:	IF IHI=IHI2 THEN GO L11;
	IHI←IHI+1;
	GO L100;

L1:	SX←SX+E11+E21;
	SY←SY+E12+E22;
	SX2←SX2+E11*E11+E21*E21;
	SY2←SY2+E12*E12+E22*E22;
	SXY←SXY+E11*E12+E21*E22;
	IF IHI1<0 THEN GO L101;
L11:	CN←2*(ABS(IHI-ILO)+1);
	F←SX*SY-CN*SXY;
	G←SY*SY-SX*SX+CN*(SX2-SY2);
	IF ABS G <500.* ABS F  THEN GO L4;
	IF ABS CY > ABS CX  THEN GO L40;
	S1←1.;
	S2←0.;
	S3←-SX/CN;
	S4←SX2+S3*SX;
	GO L51;

L40:	S1←0.;
	S2←1.;
	S3←-SY/CN;
	S4←SY2+S3*SY;
L51:	IF S4>CN*TOL THEN RETURN(RRET);
	RRET←S4/CN;
	CX←S1;
	CY←S2;
	CC←S3;
	GO L52;

L4:	G←0.5*G/F;
	ROT←SQRT(G*G+1.);
	S1←1.;
	S2←G+ROT;
	S3←G-ROT;
	SQ1←1.+S2*S2;
	SQ2←1.+S3*S3;
	G←CX+CY*S2;
	SQT←CX+CY*S3;
	IF G*G*SQ2≥SQT*SQT*SQ1 THEN GO L50;
	S2←S3;
	SQ1←SQ2;
L50:	S3←-(SX+SY*S2)/CN;
	S4←(SX2+S2*(S2*SY2+2.*SXY+S3*SY)+S3*SX)/SQ1;
	IF S4>CN*TOL THEN RETURN(RRET);
	SQT←1./SQRT(SQ1);
	CX←S1*SQT;
	CY←S2*SQT;
	CC←S3*SQT;
	RRET←S4/CN;
L52:	IF ABS CC <0.001 THEN CC←0.;
	RETURN(RRET);
	END "SUMS";
_ MALIN, ANGDIR, ANGLE;

_	Stores various line-characteristics in the initial fit.;

INTERNAL SIMPLE PROCEDURE MALIN(REAL CX,CY,CC; INTEGER LL;
			REAL SQD; INTEGER ILO,IHI);
	BEGIN "MALIN"
	LABEL L102,L104,L103;
	INTEGER IV1,IV2;
	IV2←2*LL;
	IV1←IV2-1;
	IF CX≠0 THEN GO L102;
	XLCOR[IV1]←EAX[ILO];
	YLCOR[IV1]←-CC;
	XLCOR[IV2]←EBX[IHI];
	YLCOR[IV2]←-CC;
	GO L103;

L102:	IF CY≠0 THEN GO L104;
	XLCOR[IV1]←-CC;
	YLCOR[IV1]←EAY[ILO];
	XLCOR[IV2]←-CC;
	YLCOR[IV2]←EBY[IHI];
	GO L103;

L104:	YLCOR[IV1]←CX*CX*EAY[ILO]-CY*(CX*EAX[ILO]+CC);
	XLCOR[IV1]←-(CY*YLCOR[IV1]+CC)/CX;
	YLCOR[IV2]←CX*CX*EBY[IHI]-CY*(CX*EBX[IHI]+CC);
	XLCOR[IV2]←-(CY*YLCOR[IV2]+CC)/CX;
L103:	CXL[LL]←CX;
	CYL[LL]←CY;
	CCL[LL]←CC;
	ANGLEN(LL);
	LEDG1[LL]←ILO;
	LEDG2[LL]←IHI;
	RETURN;
	END "MALIN";


_	Finds angular direction from positive abscissa to
	the vector (dx,dy).;

SIMPLE REAL PROCEDURE ANGDIR(REAL DX,DY);
	RETURN(IF (A1←57.29*AMOD(ATAN2(DY,DX)+6.2832,6.2832))≥360. THEN 0. ELSE A1);

_	Returns angle from vector (DX1,DY1) to vector (DX2,DY2);

INTERNAL SIMPLE REAL PROCEDURE ANGLE(REAL DX1,DY1,DX2,DY2);
	RETURN(AMOD(ANGDIR(DX2,DY2)-ANGDIR(DX1,DY1)+360.,360.));
_ ANGLEN, LINFIT;
_	Computes angle and length for a line.;

INTERNAL SIMPLE PROCEDURE ANGLEN(INTEGER LL);
	BEGIN "ANGLEN" INTEGER IV2;
	REAL DX,DY;
	IV2←2*LL;
	DX←XLCOR[IV2]-XLCOR[IV2-1];
	DY←YLCOR[IV2]-YLCOR[IV2-1];
	ANGARG[LL]←ANGDIR(DX,DY);
	RLEN[LL]←SQRT(DX*DX+DY*DY);
	RETURN;
	END "ANGLEN";


_	Responsible for the initial line-fit. Fits lines in as
	long successive segments as possible. Sets LCREDE ← LDATE.
	Deletes short lines, based on ILFL and RMLG. Initializes vertices.
	Returns 0 when everything is OK, otherwise the number of
	the edge-pair, at which the overflow occurred in line-space.
	Storage is assumed initialized, so start line-fit!;

INTERNAL SIMPLE INTEGER PROCEDURE LINFIT;
	BEGIN "LINFIT"
	LABEL L1000,L1010,L1003,L1001,L1,L10,L1002,L100,L101,L11,L110,
		L2,L3,L30,L4,L2001;
	INTEGER IV,ILO,IP,IHI,IRET,IRI,IW,ID,IL,IH,INOL,IV1,IV2,IDUM;
	REAL SQD,SQDD,SHR,SHR1,X,Y,CX,CY,CC,SHRIN;

	ILO←1;
	IP←1;
	IHI←1;
	NOL←1;
	IRET←0;
L1000:	IF ILO>NOEPA THEN GO L2001;
	IF NOL≤NOBAL THEN GO L1010;
	NOBAL←30+NOEPA*NOBAL/ILO;
	RETURN(ILO);

L1010:	IRI←1;
	SQD←SUMS(ILO,IHI,CX,CY,CC,RMSD);
L1003:	IW←IHI;
L1001:	CASE LE[IW] OF BEGIN ; GO L1; GO L2; GO L3; GO L4 END;
L1:	IF IRI>0 THEN GO L11;
L10:	ID←IHI-ILO+1;
	IF ID≥ILLL THEN GO L100;
	ILO←IP+1;
L1002:	IHI←ILO;
	IP←ILO;
	GO L1000;

L100:	MALIN(CX,CY,CC,NOL,SQD,ILO,IHI);
_ LINFIT cont;

_	Last line is stored. See if it can be fused with the previous
	line (iteratively).;

	IF NOL<2∨¬ICON(LEDG2[NOL-1],ILO,1) THEN GO L101;

_	Yes, fuse last two lines, provided mean square deviation is OK!;

	ILO←LEDG1[NOL-1];
	SQD←SUMS(ILO,-IHI,CX,CY,CC,RMSD);
	IF SQD=-100. THEN GO L101;
	NOL←NOL-1;
	GO L100;

L101:	ILO←IHI+1;
	NOL←NOL+1;
	GO L1002;

_ 	Change direction and pick up edges there, if possible.;

L11:	IRI←-1;
L110:	IF ILO=1∨ILO≥IHI∨LE[ILO] LAND 1∨NOL≠1∧LEDG2[NOL-1]=ILO-1 THEN GO L10;
	IL←ILO-1;
	IF DNEW(EAX[IL],EAY[IL],EBX[IL],EBY[IL],CX,CY,CC)>RDNP THEN GO L10;
	SQDD←SUMS(IHI,IL,CX,CY,CC,RMSD);
	IF SQDD=-100. THEN GO L10;
	ILO←IL;
	IW←ILO;
	SQD←SQDD;
	GO L1001;

L2:	IF IRI≤0 THEN GO L110 ELSE GO L11;
L3:	IF IRI≤0 THEN GO L10;
L30:	IF IHI=NOEPA THEN GO L11;
	IH←IHI+1;
	IF DNEW(EAX[IH],EAY[IH],EBX[IH],EBY[IH],CX,CY,CC)>RDNP THEN GO L11;
	SQDD←SUMS(ILO,IH,CX,CY,CC,RMSD);
	IF SQDD=-100. THEN GO L11;
	IHI←IH;
	SQD←SQDD;
	GO L1003;

L4:	IF IRI≤0 THEN GO L110 ELSE GO L30;
L2001:	IFREEL←NOL;
	NOL←NOL-1;

_	The line-fit is now completed. Date the lines (in LCREDE).;

	LOOP(IL,1,NOL,1) LCREDE[IL]←LDATE;
	IFREEL←NOL+1;
_ LINFIT cont, DNEW, ICON;

_	Now (!) delete short lines, based on ILFL and RMLG, and initialize
	vertices at ends of good lines, and shorten such lines by
	SHRINK at each end, to minimize overshoot.;

	INOL←NOL;
	NOV←0;
	SHRIN←1.5 MAX (SHRINK*RDEP);
	LOOP(IL,1,INOL,1)
	IF LEDG2[IL]-LEDG1[IL]+1<ILFL∨RLEN[IL]<RMLG+0.*SHRIN THEN RETLIN(IL)
		ELSE BEGIN "LPU1"

_	The line is good. Initialize its end-vertices,
	and shrink the line.;

		SHR1←SHRIN MIN (0.25*RLEN[IL]);
		SHR←SHR1/RLEN[IL];
		RLEN[IL]←RLEN[IL]-2.0*SHR1;
		IV2←2*IL;
		IV1←IV2-1;
		X←XLCOR[IV1];
		Y←YLCOR[IV1];
		XLCOR[IV1]←X+SHR*(XLCOR[IV2]-X);
		YLCOR[IV1]←Y+SHR*(YLCOR[IV2]-Y);
		XLCOR[IV2]←XLCOR[IV2]+SHR*(X-XLCOR[IV2]);
		YLCOR[IV2]←YLCOR[IV2]+SHR*(Y-YLCOR[IV2]);
		LOOP(IV,IV1,IV2,1) IDUM←MSCVCO(IV,0,1);
	        END "LPU1";
	RETURN(IRET);
	END "LINFIT";

_	Computes mean distance from projected line to new point-pair.;

INTERNAL SIMPLE REAL PROCEDURE DNEW(REAL X1,Y1,X2,Y2,CX,CY,CC);
	RETURN(0.5*(ABS(X1*CX+Y1*CY+CC)+ABS(X2*CX+Y2*CY+CC)));


_	Finds out whether two edge-pairs are in the same connected set,
	and at a distance ≤ IDIS (i.e. number of edge-pairs in between).;

INTERNAL SIMPLE INTEGER PROCEDURE ICON(INTEGER L12,L21,IDIS);
	BEGIN "ICON"
	INTEGER I3;
	IF L21-L12-1>IDIS∨LE[L12]<3∨LE[L21]≠2∧LE[L21]≠4 THEN RETURN(0);
	IF L21-L12>1 THEN LOOP(I3,L12+1,L21-1,1) IF LE[I3]≠4 THEN RETURN(0);
	RETURN(1);
	END "ICON";
_ SORTED;

_	This program sorts edge-pairs into connected subsets.;

INTERNAL PROCEDURE SORTED;
	BEGIN "SORTED"
	INTEGER IW,I8,NEXT,IP,IWEAK;
	REAL GRAV,RDEP2,DDQ,FAK,A1,A2,X,Y,DXO,DYO,DXN,DYN,DX1,DY1,DXY1,
		DXY2,D2,WEAK;
	SAFEX REAL ARRAY FAX,FAY,FBX,FBY[1:NOEPA];
	SAFEX INTEGER ARRAY IFO,IBA[1:NOEPA];

	TELL("edge-sorting");

_	Save edge-arrays and zero LE.
	Set up center coordinates in EAX and EAY.
	EBX will later contain the IBA-distance, EBY the IFO-distance.;

	ARRBLT(FAX[1],EAX[1],NOEPA);
	ARRBLT(FAY[1],EAY[1],NOEPA);
	ARRBLT(FBX[1],EBX[1],NOEPA);
	ARRBLT(FBY[1],EBY[1],NOEPA);
	LOOP(IW,1,NOEPA,1)
		BEGIN
		LE[IW]←0;
		EAX[IW]←0.5*(FAX[IW]+FBX[IW]);
		EAY[IW]←0.5*(FAY[IW]+FBY[IW]);
		EBX[IW]←1000000.;
		EBY[IW]←1000000.;
		END;
_	Establish backward and forward links.;

	GRAV←RDEP*(2.01+RDDP);
 	RDEP2←4.*RDEP↑2;
	DDQ←RDEP2*RDDP↑2;
	FAK←RDDP↑2/((15.1/(16.0*COS(RMEDA/114.58)↑4-0.9))↑2-1.0);
	A1←RDEP2↑2/16.0;
	A2←0.9*A1;
	A1←15.1*A1;
	IW←1;
	LOOP(IW,1,NOEPA-1,1)
		BEGIN "LP1" LABEL L1;

_		Get center point and pair-vector (for speed).;	

		X←EAX[IW];
		Y←EAY[IW];
		DXO←FBX[IW]-FAX[IW];
		DYO←FBY[IW]-FAY[IW];
		LOOP(NEXT,IW+1,NOEPA,1)
			BEGIN "LP100" LABEL L100,L101,L102,L103;
_ SORTED cont.;
_			Is new pair inside window?;

			IF ABS(X-EAX[NEXT])>GRAV∨ABS(Y-EAY[NEXT])>GRAV
				THEN GO L100;

_			Yes, it is. Compute directed distances and
			update minima. First find vector for new pair.;

			DXN←FBX[NEXT]-FAX[NEXT];
			DYN←FBY[NEXT]-FAY[NEXT];
			DX1←3.*(FBX[IW]-FAX[NEXT])+FBX[NEXT]-FAX[IW];
			DY1←3.*(FBY[IW]-FAY[NEXT])+FBY[NEXT]-FAY[IW];
			DXY1←DX1↑2+DY1↑2;
			DXY2←(DX1-4.*(DXO+DXN))↑2+(DY1-4.*(DYO+DYN))↑2;

_			Directed distances too large already?;

			IF DXY1>DDQ∧DXY2>DDQ THEN GO L100;

_			No. Compute, and add, correction for
			angular difference.;

			D2←1.0 MAX (A1/(0.001 MAX (((DXO+DXN)↑2+
				(DYO+DYN)↑2)↑2-A2)));
			D2←RDEP2*FAK*(D2↑2-1.0);
			DXY1←DXY1+D2;
			DXY2←DXY2+D2;


_			Now go through minimum values, and update if
			necessary.;

			IF DXY1>DDQ THEN GO L101;
			IF DXY1>EBY[IW] THEN GO L102;

_			New minimum for old forward.;

			IFO[IW]←NEXT;
			EBY[IW]←DXY1;
L102:			IF DXY1>EBX[NEXT] THEN GO L101;

_			New minimum for new backward.;

			IBA[NEXT]←IW;
			EBX[NEXT]←DXY1;
L101:			IF DXY2>DDQ THEN GO L100;
			IF DXY2>EBY[NEXT] THEN GO L103;
_ SORTED cont.;
_			New minimum for new forward.;

			IFO[NEXT]←IW;
			EBY[NEXT]←DXY2;
L103:			IF DXY2>EBX[IW] THEN GO L100;

_			New minimum for old backward.;

			IBA[IW]←NEXT;
			EBX[IW]←DXY2;
L100:	                END "LP100";
L1:		END "LP1";

_	At this point, all edge-pairs have been equipped with both
	backward and forward pointers (not necessarily reciprocated).;

_	Clean up the linkages, and break up loops (although very
	unlikely) at their weakest link.;

	LOOP(I8,1,NOEPA,1)
		BEGIN "LP8" LABEL L82,L84,L80,L81,L83,L8;
		IF LE[I8]≠0 THEN GO L8;
		WEAK←0.;
		IW←I8;
L82:		LE[IW]←1;
		NEXT←IFO[IW];

_		Chain continues?;

		IF NEXT=0∨IBA[NEXT]≠IW THEN GO L80;

_		Yes, step next.;

		IF EBY[IW]≤WEAK THEN GO L84;

_		New maximum for weak link.;

		IWEAK←IW;
		WEAK←EBY[IW];

L84:		IW←NEXT;

_		Do we have a loop?;

		IF IW≠I8 THEN GO L82;
_ SORTED cont.;
_		Yes, we have a loop. Break it at the weakest link.;

		IBA[IFO[IWEAK]]←0;
		IFO[IWEAK]←0;
		GO L8;

_		No, there is a break. Reverse.;

L80:		IFO[IW]←0;
		IW←I8;
L81:		NEXT←IBA[IW];

_		Chain continues?;

		IF NEXT=0∨IFO[NEXT]≠IW THEN GO L83;

_	Yes, step next.;

		IW←NEXT;
		LE[IW]←1;
		GO L81;

_	Break in the backward linkage. End of chain.;

L83:		IBA[IW]←0;
L8:	        END "LP8";

_	The following recopies arrays according to connectivity.;

	IP←1;
	LOOP(IW,1,NOEPA,1)
		BEGIN "LP5" LABEL L5,L6,L7;
		IF IBA[IW]≠0 THEN GO L5;
		NEXT←IW;
		LE[IP]←1;
L7:		EAX[IP]←FAX[NEXT];
		EAY[IP]←FAY[NEXT];
		EBX[IP]←FBX[NEXT];
		EBY[IP]←FBY[NEXT];
		NEXT←IFO[NEXT];
		IF NEXT=0 THEN GO L6;
		LE[IP]←LE[IP]+2;
		IP←IP+1;
		LE[IP]←2;
		GO L7;

L6:		IP←IP+1;
L5:		END "LP5";
	UNTELL;
	END "SORTED";
_ WEIGHV;

_	Computes the weighted least-square coordinates, (x,y), and the
	weight ( ← sum of square-roots of lengths of lines involved),
	for the c.v. ICV, also counting inactive lines iff ICV<0.
	If ICV>'7777 then ICV is interpreted as "number of lines" (LSH 24)
	in the array of line-numbers (VERLIN).;

INTERNAL SIMPLE PROCEDURE WEIGHV(INTEGER ICV; REFERENCE REAL X,Y,WE);
	BEGIN "WEIGHV"
	LABEL L1,L2;
	INTEGER ISV,LL,IND,NU;
	REAL SA2,SB2,SAB,SAC,SBC,SX,SY,W,DS;
	IND←1;
	NU←0;
	IF ICV>'7777 THEN NU←ICV LSH -24
		     ELSE IF¬(ISV←ABS LVNEXT(ICV,9)) THEN RETURN;
	SA2←SB2←SAB←SAC←SBC←SX←SY←WE←0.;

_	OK, that was the initialization part. Now do it!;

L1:	LL← IF NU THEN VERLIN[IND] ELSE (ISV+1)%2;
	W←SQRT(RLEN[LL]);
	WE←WE+W;
	IF ¬NU THEN BEGIN SX←SX+W*XLCOR[ISV]; SY←SY+W*YLCOR[ISV] END;
	SA2←SA2+W*CXL[LL]↑2;
	SB2←SB2+W*CYL[LL]↑2;
	SAB←SAB+W*CXL[LL]*CYL[LL];
	SAC←SAC+W*CXL[LL]*CCL[LL];
	SBC←SBC+W*CYL[LL]*CCL[LL];
	IF NU THEN IF (IND←IND+1)≤NU THEN GO L1 ELSE
	      ELSE BEGIN ISV←ABS LVNEXT(0,9); IF ISV>0 THEN GO L1 END;

_	The various sums are now computed. Find coordinates.;

	DS←SAB*SAB-SA2*SB2;

_	DS is a measure of the covariance of x- and y-coordinates, so that
	DS tends to be small for collinear lines. Hence the special case.;

	IF ABS DS > 0.1 THEN GO L2;
	X←IF NU THEN -1.0 ELSE SX/WE;
	Y←SY/WE;
	RETURN;

L2:	Y←(SA2*SBC-SAC*SAB)/DS;
	X←-(SAC+Y*SAB)/SA2;
	RETURN;
	END "WEIGHV";
_ MALI, PLDIS;

_	Finds equation and other information for the inserted line LL.;

INTERNAL SIMPLE PROCEDURE MALI(INTEGER LL; REAL X1,Y1,X2,Y2);
	BEGIN "MALI"
	INTEGER IV2;
	IV2←2*LL;
	XLCOR[IV2-1]←X1;
	YLCOR[IV2-1]←Y1;
	XLCOR[IV2]←X2;
	YLCOR[IV2]←Y2;
	LEKV(X1,Y1,X2,Y2,CXL[LL],CYL[LL],CCL[LL]);
	ANGLEN(LL);
	RETURN;
	END "MALI";


_	Finds the shortest squared distance, R, from point (X,Y) to
	line I, and the corresponding coordinates, (XL,YL), on the
	line. IW ← 1 (else 0) iff (XL,YL) is outside the line segment.
	This routine is used in the insertion package. Assumes the
	topological connectivity as reflected in the line-coordinates.;

INTERNAL SIMPLE PROCEDURE PLDIS(REAL X,Y; INTEGER I;
		REFERENCE REAL XL,YL,R; REFERENCE INTEGER IW);
	BEGIN "PLDIS"
	LABEL L2;
	INTEGER IV;
	REAL AK;
	IW←0;
	IV←2*I-1;
	AK←1000.;
	IF CYL[I]≠0 THEN AK←-CXL[I]/CYL[I];
	YL←(YLCOR[IV]+AK*(AK*Y-XLCOR[IV]+X))/(1.+AK*AK);
	XL←X-AK*(YL-Y);
	R←(X-XL)↑2+(Y-YL)↑2;
	IF ABS AK > 1. THEN GO L2;
	IF (XL-XLCOR[IV])*(XL-XLCOR[IV+1])≥-1. THEN IW←1;
	RETURN;

L2:	IF (YL-YLCOR[IV])*(YL-YLCOR[IV+1])≥-1. THEN IW←1;
	RETURN;
	END "PLDIS";
_ LDIST, SKAR1, LEKV;

_	Measures distance (signed) from (X,Y) to line L.;

INTERNAL SIMPLE REAL PROCEDURE LDIST(REAL X,Y; INTEGER L);
	RETURN(X*CXL[L]+Y*CYL[L]+CCL[L]);


_	Finds intersection of (X1,Y1)-(X2,Y2) with line LL, after 
	(X1,Y1) on the ray only. Returns point of intersection, (X,Y),
	and distance from (X1,Y1), RSQ.
	RSQ ← 900000. iff there is no such intersection.;

INTERNAL SIMPLE PROCEDURE SKAR1(REAL X1,Y1,X2,Y2; INTEGER LL; 
	        REFERENCE REAL X,Y,RSQ);
	BEGIN "SKAR1"
	INTEGER IV,IDUM,I1,I2,I3,I4;
	REAL R2;
	IV←2*LL;
	IDUM←KARN(X1,Y1,X2,Y2,XLCOR[IV-1],YLCOR[IV-1],
		  XLCOR[IV],YLCOR[IV],X,Y,I1,I2,I3,I4,RSQ,R2,1,0.);
	RSQ←(X1-X)↑2+(Y1-Y)↑2;
	IF I3=0∨I3=1∨I3=-1∧RSQ<0.25 THEN RSQ←900000.;
	RETURN;
	END "SKAR1";


_	Finds the normalized line-equation.;

INTERNAL SIMPLE PROCEDURE LEKV(REAL X1,Y1,X2,Y2; REFERENCE REAL A,B,C);
	BEGIN "LEKV"
	REAL R;
	IF ABS(Y1-Y2)≤0.01 THEN BEGIN A←0.;B←1.;C←-0.5*(Y1+Y2);RETURN END;
	A←1.;
	B←(X2-X1)/(Y1-Y2);
	R←SQRT(A↑2+B↑2);
	IF ABS (A←A/R) < 0.002  THEN A←0.;
	IF ABS (B←B/R) < 0.002 THEN B←0.;
	IF ABS (C←-A*X1-B*Y1) < 0.001 THEN C←0.
	END "LEKV";
_ REKOP, INREK;

_	Sets up transformation data from internal representation into
	a rectangular operator with the line (X1,Y1)-(X2,Y2) as axis,
	and of width RR. RL returns ratio of length of rectangle to
	width of rectangle. The origin of this new coordinate system
	is positioned at the centre of gravity of the rectangle.
	The new X-axis is directed to (X2,Y2), and the Y-axis goes
	out perpendicularly on the left.;

INTERNAL SIMPLE PROCEDURE REKOP(REAL X1,Y1,X2,Y2,RR; REFERENCE REAL RL);
	BEGIN "REKOP"
	REAL DX2,DY2,DX3,DY3,Q2,Q3,X,Y,X3,Y3;
	X←0.5*(X1+X2);
	Y←0.5*(Y1+Y2);
	RL←SQRT((X1-X2)↑2+(Y1-Y2)↑2)/RR;
	X3 ← X-(Y2-Y)/RL;
	Y3 ← Y+(X2-X)/RL;
	DY2←SIGN(ABS(DX2←Y2-Y) MAX 0.000001,DX2);
	DY3←SIGN(ABS(DX3←Y3-Y) MAX 0.000001,DX3);
	DX2←X2-X;
	DX3←X3-X;
	Q2←DX2/DY2;
	Q3←DX3/DY3;
	A11←1./(DX2-Q3*DY2);
	A21←1./(DX3-Q2*DY3);

_	Note that e.g. DX2-Q3*DY2=0 would imply parallelity of the axes.
	Also e.g. DX2 and DX3 cannot simultaneously be 0.;

	IF ABS(Q3)<.000001 THEN Q3←IF Q3≥0 THEN .000001 ELSE -.000001;
	IF ABS(Q2)<.000001 THEN Q2←IF Q2≥0 THEN .000001 ELSE -.000001;
	A12←1./(DY2-DX2/Q3);
	A22←1./(DY3-DX3/Q2);
	X00←X;
	Y00←Y
	END "REKOP";


_	Returns T (else F) iff (X,Y) is inside current rectangular operator.;

INTERNAL SIMPLE INTEGER PROCEDURE INREK(REAL X,Y);
	RETURN(ABS(A11*(X-X00)+A12*(Y-Y00))≤1.∧
		ABS(A21*(X-X00)+A22*(Y-Y00))≤1.);
END "LINE";